Brms_growth_lobe_number_Gompertz_4p

Go to the end of this file for the parameter estimates of the best fitted model on phenotype lobe number

fitting non-linear model with BRMS (Bayesian Regression Models using Stan)

lobe.n <- read_csv("lobe.n.csv", col_types = list(col_character(), col_date(), col_double(), col_double()))
lobe.n 

sample.list <- unique(lobe.n$ID)

166 individuals at 4 time points

ggplot(data=lobe.n, aes(x=days, y=lobe_n, group=ID)) +
    geom_line(alpha=.1) + 
    geom_point(size=.1, alpha=.05) +
    ggtitle("lobe_number by days") 

#### Try Gompertz (4-parameter) Model (see http://www.pisces-conservation.com/growthhelp/index.html?weibul.html) by using BRMS

First set up the formula

gompertz_4p.bf1 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k*(days - I))), 
  Hmax + Hmin + k + I ~ 1, #general model, paramters do not vary for individuals
  nl=TRUE)

Hmax, asymptotic height at which growth is zero
Hmin, lower asymptotic height
k, growth rate
I, time at the inflection point

stat.data <- function(data) {
  data %>% 
  group_by(days) %>%
  summarize(meadian=median(lobe_n),
            max=max(lobe_n),
            min=min(lobe_n),
            sd=sd(lobe_n))
}
stat.data(lobe.n)

Priors. Hmin and Hmax use median at start and end dates.

prior1 <- c(prior(normal(11,4), nlpar="Hmax"), 
            prior(normal(0,3), nlpar="Hmin"),
            prior(normal(1,1), nlpar="k"), 
            prior(normal(150,10), nlpar="I"))
fit1 <- brm(formula=gompertz_4p.bf1,
            data=lobe.n,
            prior=prior1)

summary(fit1, waic=TRUE, R2=TRUE)

The model has not converged (some Rhats are > 1.1). Do not analyse the results! 
We recommend running more iterations and/or setting stronger priors.There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I))) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         I ~ 1
   Data: lobe.n (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 3598.29; R2 = 0.56
 
Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    10.99      0.46    10.15    11.87          4 1.35
Hmin_Intercept     1.08      0.29     0.53     1.65        122 1.04
k_Intercept        0.29      0.35     0.06     0.93          2 8.46
I_Intercept      136.71      1.53   134.93   140.09          5 1.27

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.63      0.10     3.45     3.82        106 1.03

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit1)

pairs(fit1)

keep k positive. Or more tightly constrain the priors on Hmax and Hmin.

summary(fit2, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I))) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         I ~ 1
   Data: lobe.n (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 3595.49; R2 = 0.56
 
Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    11.17      0.39    10.44    12.01       1218 1.00
Hmin_Intercept     1.07      0.29     0.49     1.62       2396 1.00
k_Intercept        0.09      0.04     0.06     0.21        464 1.01
I_Intercept      137.24      1.37   134.88   140.09       1536 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.61      0.10     3.42     3.82       3255 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit2)

pairs(fit2)

to make k a little bit reasonable

gompertz_4p.bf2 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10)*(days - I))), 
  Hmax + Hmin + k + I ~ 1, #general model, paramters do not vary for individuals
  nl=TRUE)

summary(fit3, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         I ~ 1
   Data: lobe.n (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 3595.29; R2 = 0.56
 
Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    11.13      0.37    10.43    11.89       1847 1.00
Hmin_Intercept     1.01      0.30     0.42     1.58       2465 1.00
k_Intercept        0.84      0.23     0.56     1.43       1332 1.00
I_Intercept      136.80      1.30   134.40   139.49       2425 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.61      0.10     3.42     3.81       3699 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit3)

pairs(fit3)

What does the fit look like?

newdata <- data.frame(days=seq(min(lobe.n$days), max(lobe.n$days),1))
fit3.fitted <-  cbind(newdata, fitted(fit3, newdata)) %>% as.tibble() %>%
  rename(lobe_n=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')

plot

pl <- ggplot(aes(x=days, y=lobe_n),data=NULL)
pl <- pl + geom_line(aes(group=ID), alpha=.1, data=lobe.n)
pl + geom_line(color="skyblue", lwd=1.5, data=fit3.fitted)

Now try adding random effects for model parameters:

What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:

gompertz_4p.fn <- function(Hmax, Hmin, k, I, days) {
    Hmin + (Hmax - Hmin) * exp(-exp(-(k/10)*(days - I))) 
}
for(I1 in seq(100,170,10)) {
  for(k1 in seq(0,1,.2)) {
  tmp.lobe.n <- gompertz_4p.fn(Hmax=11,
                               Hmin=0,
                               k=k1,
                               I=I1,
                              days=newdata$days)
  abc <- data.frame(newdata$days, tmp.lobe.n)
  p <- ggplot(data=abc, aes(newdata$days, tmp.lobe.n)) +
    geom_line() + ylim(0,20) + ggtitle(paste0("I=",I1," k=",k1))
  print(p)
  }
}

First try with only fixing Hmin

gompertz_4p.bf3 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10)*(days - I))), 
  Hmax + k + I ~ (1|ID), # vary for individuals
  Hmin ~ 1, # do not vary per individual
  nl=TRUE) 

summary(fit5, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))) 
         Hmax ~ (1 | ID)
         k ~ (1 | ID)
         I ~ (1 | ID)
         Hmin ~ 1
   Data: lobe.n (Number of observations: 664) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = 2871.71; R2 = 0.91
 
Group-Level Effects: 
~ID (Number of levels: 166) 
                   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept)     2.27      0.18     1.94     2.64       1743 1.00
sd(k_Intercept)        0.97      0.20     0.58     1.33        557 1.01
sd(I_Intercept)       17.22      1.14    15.14    19.63       1309 1.00

Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    11.31      0.23    10.85    11.76        915 1.01
k_Intercept        2.24      0.23     1.83     2.71       4562 1.00
I_Intercept      138.75      1.52   135.79   141.77       1390 1.00
Hmin_Intercept     0.42      0.12     0.19     0.66      10000 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.65      0.07     1.52     1.80       1493 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

get fitted values

b <- lobe.n %>% ungroup()
fit5.fitted <- cbind(b, fitted(fit5)) %>% as.tibble()
fit5.fitted

plot

plot.fitted(fit5.fitted)

plot fitted vs actual:

plot.fitted.actual <- function(fn) {
  
  r_squared <- summary(lm(Estimate ~ lobe_n, data=fn))$adj.r.squared
  r_squared <- round(r_squared, digits = 4)

  fn %>%
    mutate(days=as.factor(days)) %>%
    ggplot(aes(x=lobe_n, y=Estimate, shape=days, color=days)) +
    geom_point() +
    geom_abline(intercept = 0, slope=1) +
    scale_x_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20)) +
    ggtitle(paste0("R2 = ",r_squared))

}

total_sum_of_squared_residuals <- function(fn) {anova(lm(Estimate ~ lobe_n, data=fn))[2,2]}

plot.fitted.actual(fit5.fitted)

cv

SSR.fit <- total_sum_of_squared_residuals(fit5.fitted)
waic.fit <- round(waic(fit5)$waic, digits=2)
kfoldic.fit <- round(kfold(fit5)$kfoldic, digits=2)

lobe.n.mean

lobe.n.mean <- read_csv("lobe.n.mean.csv")
lobe.n.mean

 lapply(1:6, function(i) {
       ggplot(aes(x=days, group=ID), data=lobe.n.mean) +
         geom_line(aes(y=lobe_n), color="red") +   #red: modified data
         geom_line(aes(y=lobe_n_raw), color="black") +
          facet_wrap_paginate(~ID, ncol = 6, nrow = 5,  page =i)
}
)

Try Gompertz with 4-parameter (see http://www.pisces-conservation.com/growthhelp/index.html?weibul.html) by using BRMS

First set up the formula

gompertz_4p.bf1 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k*(days - I))), 
  Hmax + Hmin + k + I ~ 1, #general model, paramters do not vary for individuals
  nl=TRUE)

stat.data(lobe.n.mean)

stat.data(lobe.n)

Priors. Hmin and Hmax use median at start and end dates.

summary(fit1.mean, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I))) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         I ~ 1
   Data: lobe.n.mean (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 3594.58; R2 = 0.6
 
Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    12.58      0.46    11.77    13.55       1520 1.00
Hmin_Intercept     0.84      0.41    -0.13     1.51       1310 1.00
k_Intercept        0.06      0.01     0.04     0.08       1750 1.00
I_Intercept      136.44      1.62   133.03   139.39       1741 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.61      0.10     3.43     3.82       2901 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit1.mean)

pairs(fit1.mean)

summary(fit2.mean, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I))) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         I ~ 1
   Data: lobe.n.mean (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 3594.93; R2 = 0.6
 
Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    12.58      0.49    11.75    13.68       1262 1.00
Hmin_Intercept     0.83      0.46    -0.24     1.52        977 1.01
k_Intercept        0.06      0.01     0.04     0.08       1320 1.00
I_Intercept      136.43      1.78   132.72   139.66       1325 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.61      0.10     3.43     3.81       3408 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit2.mean)

pairs(fit2.mean)

gompertz_4p.bf2 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))),
  Hmax + Hmin + k + I ~ 1, 
  nl=TRUE)

summary(fit3.mean, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         I ~ 1
   Data: lobe.n.mean (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 3594.99; R2 = 0.6
 
Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    12.58      0.48    11.76    13.68       1377 1.00
Hmin_Intercept     0.81      0.47    -0.39     1.52        767 1.00
k_Intercept        0.60      0.11     0.41     0.83       1369 1.00
I_Intercept      136.34      1.82   131.76   139.45       1032 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.61      0.10     3.43     3.81       2695 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit3.mean)

pairs(fit3.mean)

What does the fit look like?

newdata <- data.frame(days=seq(min(lobe.n$days), max(lobe.n$days),1))
fit3.mean.fitted <-  cbind(newdata, fitted(fit3.mean, newdata)) %>% as.tibble() %>%
  rename(lobe_n=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')

plot

ggplot(aes(x=days, y=lobe_n),data=NULL) +
  geom_line(aes(group=ID), alpha=.1, data=lobe.n) +
  geom_line(color="skyblue", lwd=1.5, data=fit3.mean.fitted)

What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:

lobe.n.fn <- function (Hmax, Hmin, k, delta, days) {
  Hmax - (Hmax - Hmin) * exp(-(k/10^6) * (days^delta)) 
}
for(delta1 in seq(2,4,.25)) {
  tmp.lobe.n <- lobe.n.fn(Hmax=8, 
                          Hmin=-2, 
                          k=.29, 
                          delta=delta1, 
                          days=newdata$days)
   abc <- data.frame(newdata$days, tmp.lobe.n)
   p <- ggplot(data=abc, aes(newdata$days, tmp.lobe.n)) +
    geom_line() + ylim(0,40) + ggtitle("delta=",delta1)
   print(p)
}

for(k1 in seq(0,1,.25)) {
  tmp.lobe.n <- lobe.n.fn(Hmax=8, 
                          Hmin=-2, 
                          k=k1, 
                          delta=3, 
                          days=newdata$days)
   abc <- data.frame(newdata$days, tmp.lobe.n)
   p <- ggplot(data=abc, aes(newdata$days, tmp.lobe.n)) +
    geom_line() + ylim(0,40) + ggtitle("k=",k1)
  print(p)
}

gompertz_4p.bf3 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))),
  Hmax + k + I ~ (1|ID), # vary for individuals
  Hmin ~ 1, # do not vary per individual
  nl=TRUE)

summary(fit5.mean, waic=TRUE, R2=TRUE)

There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))) 
         Hmax ~ (1 | ID)
         k ~ (1 | ID)
         I ~ (1 | ID)
         Hmin ~ 1
   Data: lobe.n.mean (Number of observations: 664) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = 1803.34; R2 = 0.99
 
Group-Level Effects: 
~ID (Number of levels: 166) 
                   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept)     2.88      0.17     2.56     3.25       1208 1.01
sd(k_Intercept)        0.94      0.10     0.75     1.16       1959 1.00
sd(I_Intercept)       18.12      1.10    16.06    20.37       1284 1.01

Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    12.63      0.24    12.14    13.10        729 1.01
k_Intercept        1.57      0.12     1.35     1.83       2552 1.00
I_Intercept      136.87      1.43   134.11   139.73        465 1.01
Hmin_Intercept     0.06      0.05    -0.04     0.16      10000 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.67      0.03     0.61     0.72       2403 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

get fitted values

b <- lobe.n.mean %>% ungroup()
fit5.mean.fitted <- cbind(b, fitted(fit5.mean)) %>% as.tibble()
fit5.mean.fitted

plot

plot.fitted(fit5.mean.fitted)

plot fitted vs actual(modified)

plot.fitted.actual(fit5.mean.fitted)

plot fitted vs actual(raw)

plot.fitted.actual.raw(fit5.mean.fitted)

SSR on raw data set

total_sum_of_squared_residuals_raw <- function(fn) {anova(lm(Estimate ~ lobe_n_raw, data=fn))[2,2]}

SSR.fit.mean.raw <- total_sum_of_squared_residuals_raw(fit5.mean.fitted)

cv

SSR.fit.mean <- total_sum_of_squared_residuals(fit5.mean.fitted)
WAIC.fit.mean <- round(waic(fit5.mean)$waic, digits=2)
KFOLDIC.fit.mean <- round(kfold(fit5.mean)$kfoldic, digits=2)

lobe.n.pmm

lobe.n.pmm <- read_csv("lobe.n.pmm.csv")
lobe.n.pmm

 lapply(1:6, function(i) {
       ggplot(aes(x=days, group=ID), data=lobe.n.pmm) +
         geom_line(aes(y=lobe_n), color="red") +   #red: modified data
         geom_line(aes(y=lobe_n_raw), color="black") +
          facet_wrap_paginate(~ID, ncol = 6, nrow = 5,  page =i)
}
)

stat.data(lobe.n.pmm)

stat.data(lobe.n.mean)

stat.data(lobe.n)

gompertz_4p.bf3 <- bf(
  lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))),
  Hmax + k + I ~ (1|ID), # vary for individuals
  Hmin ~ 1, # do not vary per individual
  nl=TRUE)

summary(fit5.pmm, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))) 
         Hmax ~ (1 | ID)
         k ~ (1 | ID)
         I ~ (1 | ID)
         Hmin ~ 1
   Data: lobe.n.pmm (Number of observations: 664) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = 2222.45; R2 = 0.97
 
Group-Level Effects: 
~ID (Number of levels: 166) 
                   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept)     2.78      0.18     2.46     3.15       1954 1.00
sd(k_Intercept)        0.90      0.12     0.69     1.17       2514 1.00
sd(I_Intercept)       17.50      1.07    15.48    19.75       1565 1.01

Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept    12.50      0.24    12.03    12.97       1614 1.00
k_Intercept        1.80      0.16     1.51     2.14       5011 1.00
I_Intercept      138.42      1.45   135.58   141.29       1015 1.00
Hmin_Intercept     0.22      0.07     0.08     0.36      10000 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.97      0.04     0.89     1.06       2727 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

get fitted values

b <- lobe.n.pmm %>% ungroup()
fit5.pmm.fitted <- cbind(b, fitted(fit5.pmm)) %>% as.tibble()
fit5.pmm.fitted

plot

plot.fitted(fit5.pmm.fitted)

plot fitted vs actual(modified):

plot.fitted.actual(fit5.pmm.fitted)

plot fitted vs actual(raw):

plot.fitted.actual.raw(fit5.pmm.fitted)

SSR on raw data set

SSR.fit.pmm.raw <- total_sum_of_squared_residuals_raw(fit5.pmm.fitted)

cv

SSR.fit.pmm <- total_sum_of_squared_residuals(fit5.pmm.fitted)
WAIC.fit.pmm <- round(waic(fit5.pmm)$waic, digits=2)
KFOLDIC.fit.pmm <- round(kfold(fit5.pmm)$kfoldic, digits=2)

which fitted model is the better?

create an empty list for models

0. SSR on raw data set

cv.list[2,1] <- SSR.fit.mean.raw
cv.list[3,1] <- SSR.fit.pmm.raw

1. SSR

cv.list[1,2] <- SSR.fit
cv.list[2,2] <- SSR.fit.mean
cv.list[3,2] <- SSR.fit.pmm

2. CV

WAIC(Widely Applicable Information Criterion) an extension of the Akaike Information Criterion (AIC) that is more fully Bayesian.

K-fold CV Data are randomly partitioned into K subsets of equal size. Then the model is refit 10 times(default), each time leaving out one of the 10 subsets.

cv.list[1,3] <- waic.fit
cv.list[2,3] <- WAIC.fit.mean
cv.list[3,3] <- WAIC.fit.pmm

cv.list[1,4] <- kfoldic.fit
cv.list[2,4] <- KFOLDIC.fit.mean
cv.list[3,4] <- KFOLDIC.fit.pmm

cv.list

parameter estimates of best fitted model on phenotype lobe number : take fit5.mean

read 
1) population (fixed) & group level (random) effects
2) parameters: Hmax, k & I of each genotype

lobe.n.best.fitted <- fit5.mean$fit
dimnames(lobe.n.best.fitted)

$iterations
NULL

$chains
[1] "chain:1" "chain:2" "chain:3" "chain:4"

$parameters
  [1] "b_Hmax_Intercept"             "b_k_Intercept"                "b_I_Intercept"               
  [4] "b_Hmin_Intercept"             "sd_ID__Hmax_Intercept"        "sd_ID__k_Intercept"          
  [7] "sd_ID__I_Intercept"           "sigma"                        "r_ID__Hmax[ID_1,Intercept]"  
 [10] "r_ID__Hmax[ID_10,Intercept]"  "r_ID__Hmax[ID_100,Intercept]" "r_ID__Hmax[ID_101,Intercept]"
 [13] "r_ID__Hmax[ID_102,Intercept]" "r_ID__Hmax[ID_103,Intercept]" "r_ID__Hmax[ID_104,Intercept]"
 [16] "r_ID__Hmax[ID_105,Intercept]" "r_ID__Hmax[ID_106,Intercept]" "r_ID__Hmax[ID_107,Intercept]"
 [19] "r_ID__Hmax[ID_109,Intercept]" "r_ID__Hmax[ID_111,Intercept]" "r_ID__Hmax[ID_112,Intercept]"
 [22] "r_ID__Hmax[ID_113,Intercept]" "r_ID__Hmax[ID_114,Intercept]" "r_ID__Hmax[ID_116,Intercept]"
 [25] "r_ID__Hmax[ID_117,Intercept]" "r_ID__Hmax[ID_118,Intercept]" "r_ID__Hmax[ID_119,Intercept]"
 [28] "r_ID__Hmax[ID_12,Intercept]"  "r_ID__Hmax[ID_120,Intercept]" "r_ID__Hmax[ID_121,Intercept]"
 [31] "r_ID__Hmax[ID_122,Intercept]" "r_ID__Hmax[ID_123,Intercept]" "r_ID__Hmax[ID_124,Intercept]"
 [34] "r_ID__Hmax[ID_125,Intercept]" "r_ID__Hmax[ID_126,Intercept]" "r_ID__Hmax[ID_127,Intercept]"
 [37] "r_ID__Hmax[ID_128,Intercept]" "r_ID__Hmax[ID_129,Intercept]" "r_ID__Hmax[ID_13,Intercept]" 
 [40] "r_ID__Hmax[ID_130,Intercept]" "r_ID__Hmax[ID_132,Intercept]" "r_ID__Hmax[ID_133,Intercept]"
 [43] "r_ID__Hmax[ID_134,Intercept]" "r_ID__Hmax[ID_135,Intercept]" "r_ID__Hmax[ID_136,Intercept]"
 [46] "r_ID__Hmax[ID_137,Intercept]" "r_ID__Hmax[ID_138,Intercept]" "r_ID__Hmax[ID_139,Intercept]"
 [49] "r_ID__Hmax[ID_14,Intercept]"  "r_ID__Hmax[ID_142,Intercept]" "r_ID__Hmax[ID_143,Intercept]"
 [52] "r_ID__Hmax[ID_144,Intercept]" "r_ID__Hmax[ID_147,Intercept]" "r_ID__Hmax[ID_148,Intercept]"
 [55] "r_ID__Hmax[ID_15,Intercept]"  "r_ID__Hmax[ID_152,Intercept]" "r_ID__Hmax[ID_153,Intercept]"
 [58] "r_ID__Hmax[ID_154,Intercept]" "r_ID__Hmax[ID_155,Intercept]" "r_ID__Hmax[ID_156,Intercept]"
 [61] "r_ID__Hmax[ID_157,Intercept]" "r_ID__Hmax[ID_158,Intercept]" "r_ID__Hmax[ID_16,Intercept]" 
 [64] "r_ID__Hmax[ID_160,Intercept]" "r_ID__Hmax[ID_161,Intercept]" "r_ID__Hmax[ID_163,Intercept]"
 [67] "r_ID__Hmax[ID_165,Intercept]" "r_ID__Hmax[ID_166,Intercept]" "r_ID__Hmax[ID_169,Intercept]"
 [70] "r_ID__Hmax[ID_17,Intercept]"  "r_ID__Hmax[ID_170,Intercept]" "r_ID__Hmax[ID_171,Intercept]"
 [73] "r_ID__Hmax[ID_172,Intercept]" "r_ID__Hmax[ID_174,Intercept]" "r_ID__Hmax[ID_175,Intercept]"
 [76] "r_ID__Hmax[ID_176,Intercept]" "r_ID__Hmax[ID_177,Intercept]" "r_ID__Hmax[ID_18,Intercept]" 
 [79] "r_ID__Hmax[ID_180,Intercept]" "r_ID__Hmax[ID_181,Intercept]" "r_ID__Hmax[ID_186,Intercept]"
 [82] "r_ID__Hmax[ID_189,Intercept]" "r_ID__Hmax[ID_19,Intercept]"  "r_ID__Hmax[ID_191,Intercept]"
 [85] "r_ID__Hmax[ID_192,Intercept]" "r_ID__Hmax[ID_197,Intercept]" "r_ID__Hmax[ID_198,Intercept]"
 [88] "r_ID__Hmax[ID_199,Intercept]" "r_ID__Hmax[ID_2,Intercept]"   "r_ID__Hmax[ID_20,Intercept]" 
 [91] "r_ID__Hmax[ID_201,Intercept]" "r_ID__Hmax[ID_202,Intercept]" "r_ID__Hmax[ID_207,Intercept]"
 [94] "r_ID__Hmax[ID_208,Intercept]" "r_ID__Hmax[ID_21,Intercept]"  "r_ID__Hmax[ID_212,Intercept]"
 [97] "r_ID__Hmax[ID_214,Intercept]" "r_ID__Hmax[ID_217,Intercept]" "r_ID__Hmax[ID_223,Intercept]"
[100] "r_ID__Hmax[ID_224,Intercept]" "r_ID__Hmax[ID_227,Intercept]" "r_ID__Hmax[ID_23,Intercept]" 
[103] "r_ID__Hmax[ID_231,Intercept]" "r_ID__Hmax[ID_234,Intercept]" "r_ID__Hmax[ID_24,Intercept]" 
[106] "r_ID__Hmax[ID_25,Intercept]"  "r_ID__Hmax[ID_26,Intercept]"  "r_ID__Hmax[ID_27,Intercept]" 
[109] "r_ID__Hmax[ID_28,Intercept]"  "r_ID__Hmax[ID_29,Intercept]"  "r_ID__Hmax[ID_3,Intercept]"  
[112] "r_ID__Hmax[ID_30,Intercept]"  "r_ID__Hmax[ID_32,Intercept]"  "r_ID__Hmax[ID_34,Intercept]" 
[115] "r_ID__Hmax[ID_35,Intercept]"  "r_ID__Hmax[ID_36,Intercept]"  "r_ID__Hmax[ID_37,Intercept]" 
[118] "r_ID__Hmax[ID_38,Intercept]"  "r_ID__Hmax[ID_39,Intercept]"  "r_ID__Hmax[ID_4,Intercept]"  
[121] "r_ID__Hmax[ID_40,Intercept]"  "r_ID__Hmax[ID_41,Intercept]"  "r_ID__Hmax[ID_42,Intercept]" 
[124] "r_ID__Hmax[ID_45,Intercept]"  "r_ID__Hmax[ID_47,Intercept]"  "r_ID__Hmax[ID_48,Intercept]" 
[127] "r_ID__Hmax[ID_49,Intercept]"  "r_ID__Hmax[ID_50,Intercept]"  "r_ID__Hmax[ID_51,Intercept]" 
[130] "r_ID__Hmax[ID_52,Intercept]"  "r_ID__Hmax[ID_53,Intercept]"  "r_ID__Hmax[ID_55,Intercept]" 
[133] "r_ID__Hmax[ID_56,Intercept]"  "r_ID__Hmax[ID_57,Intercept]"  "r_ID__Hmax[ID_58,Intercept]" 
[136] "r_ID__Hmax[ID_59,Intercept]"  "r_ID__Hmax[ID_6,Intercept]"   "r_ID__Hmax[ID_60,Intercept]" 
[139] "r_ID__Hmax[ID_61,Intercept]"  "r_ID__Hmax[ID_62,Intercept]"  "r_ID__Hmax[ID_63,Intercept]" 
[142] "r_ID__Hmax[ID_64,Intercept]"  "r_ID__Hmax[ID_65,Intercept]"  "r_ID__Hmax[ID_67,Intercept]" 
[145] "r_ID__Hmax[ID_68,Intercept]"  "r_ID__Hmax[ID_69,Intercept]"  "r_ID__Hmax[ID_7,Intercept]"  
[148] "r_ID__Hmax[ID_72,Intercept]"  "r_ID__Hmax[ID_73,Intercept]"  "r_ID__Hmax[ID_74,Intercept]" 
[151] "r_ID__Hmax[ID_75,Intercept]"  "r_ID__Hmax[ID_76,Intercept]"  "r_ID__Hmax[ID_77,Intercept]" 
[154] "r_ID__Hmax[ID_78,Intercept]"  "r_ID__Hmax[ID_79,Intercept]"  "r_ID__Hmax[ID_8,Intercept]"  
[157] "r_ID__Hmax[ID_80,Intercept]"  "r_ID__Hmax[ID_81,Intercept]"  "r_ID__Hmax[ID_82,Intercept]" 
[160] "r_ID__Hmax[ID_84,Intercept]"  "r_ID__Hmax[ID_85,Intercept]"  "r_ID__Hmax[ID_86,Intercept]" 
[163] "r_ID__Hmax[ID_87,Intercept]"  "r_ID__Hmax[ID_88,Intercept]"  "r_ID__Hmax[ID_89,Intercept]" 
[166] "r_ID__Hmax[ID_9,Intercept]"   "r_ID__Hmax[ID_90,Intercept]"  "r_ID__Hmax[ID_91,Intercept]" 
[169] "r_ID__Hmax[ID_92,Intercept]"  "r_ID__Hmax[ID_94,Intercept]"  "r_ID__Hmax[ID_95,Intercept]" 
[172] "r_ID__Hmax[ID_96,Intercept]"  "r_ID__Hmax[ID_98,Intercept]"  "r_ID__Hmax[ID_99,Intercept]" 
[175] "r_ID__k[ID_1,Intercept]"      "r_ID__k[ID_10,Intercept]"     "r_ID__k[ID_100,Intercept]"   
[178] "r_ID__k[ID_101,Intercept]"    "r_ID__k[ID_102,Intercept]"    "r_ID__k[ID_103,Intercept]"   
[181] "r_ID__k[ID_104,Intercept]"    "r_ID__k[ID_105,Intercept]"    "r_ID__k[ID_106,Intercept]"   
[184] "r_ID__k[ID_107,Intercept]"    "r_ID__k[ID_109,Intercept]"    "r_ID__k[ID_111,Intercept]"   
[187] "r_ID__k[ID_112,Intercept]"    "r_ID__k[ID_113,Intercept]"    "r_ID__k[ID_114,Intercept]"   
[190] "r_ID__k[ID_116,Intercept]"    "r_ID__k[ID_117,Intercept]"    "r_ID__k[ID_118,Intercept]"   
[193] "r_ID__k[ID_119,Intercept]"    "r_ID__k[ID_12,Intercept]"     "r_ID__k[ID_120,Intercept]"   
[196] "r_ID__k[ID_121,Intercept]"    "r_ID__k[ID_122,Intercept]"    "r_ID__k[ID_123,Intercept]"   
[199] "r_ID__k[ID_124,Intercept]"    "r_ID__k[ID_125,Intercept]"    "r_ID__k[ID_126,Intercept]"   
[202] "r_ID__k[ID_127,Intercept]"    "r_ID__k[ID_128,Intercept]"    "r_ID__k[ID_129,Intercept]"   
[205] "r_ID__k[ID_13,Intercept]"     "r_ID__k[ID_130,Intercept]"    "r_ID__k[ID_132,Intercept]"   
[208] "r_ID__k[ID_133,Intercept]"    "r_ID__k[ID_134,Intercept]"    "r_ID__k[ID_135,Intercept]"   
[211] "r_ID__k[ID_136,Intercept]"    "r_ID__k[ID_137,Intercept]"    "r_ID__k[ID_138,Intercept]"   
[214] "r_ID__k[ID_139,Intercept]"    "r_ID__k[ID_14,Intercept]"     "r_ID__k[ID_142,Intercept]"   
[217] "r_ID__k[ID_143,Intercept]"    "r_ID__k[ID_144,Intercept]"    "r_ID__k[ID_147,Intercept]"   
[220] "r_ID__k[ID_148,Intercept]"    "r_ID__k[ID_15,Intercept]"     "r_ID__k[ID_152,Intercept]"   
[223] "r_ID__k[ID_153,Intercept]"    "r_ID__k[ID_154,Intercept]"    "r_ID__k[ID_155,Intercept]"   
[226] "r_ID__k[ID_156,Intercept]"    "r_ID__k[ID_157,Intercept]"    "r_ID__k[ID_158,Intercept]"   
[229] "r_ID__k[ID_16,Intercept]"     "r_ID__k[ID_160,Intercept]"    "r_ID__k[ID_161,Intercept]"   
[232] "r_ID__k[ID_163,Intercept]"    "r_ID__k[ID_165,Intercept]"    "r_ID__k[ID_166,Intercept]"   
[235] "r_ID__k[ID_169,Intercept]"    "r_ID__k[ID_17,Intercept]"     "r_ID__k[ID_170,Intercept]"   
[238] "r_ID__k[ID_171,Intercept]"    "r_ID__k[ID_172,Intercept]"    "r_ID__k[ID_174,Intercept]"   
[241] "r_ID__k[ID_175,Intercept]"    "r_ID__k[ID_176,Intercept]"    "r_ID__k[ID_177,Intercept]"   
[244] "r_ID__k[ID_18,Intercept]"     "r_ID__k[ID_180,Intercept]"    "r_ID__k[ID_181,Intercept]"   
[247] "r_ID__k[ID_186,Intercept]"    "r_ID__k[ID_189,Intercept]"    "r_ID__k[ID_19,Intercept]"    
[250] "r_ID__k[ID_191,Intercept]"    "r_ID__k[ID_192,Intercept]"    "r_ID__k[ID_197,Intercept]"   
[253] "r_ID__k[ID_198,Intercept]"    "r_ID__k[ID_199,Intercept]"    "r_ID__k[ID_2,Intercept]"     
[256] "r_ID__k[ID_20,Intercept]"     "r_ID__k[ID_201,Intercept]"    "r_ID__k[ID_202,Intercept]"   
[259] "r_ID__k[ID_207,Intercept]"    "r_ID__k[ID_208,Intercept]"    "r_ID__k[ID_21,Intercept]"    
[262] "r_ID__k[ID_212,Intercept]"    "r_ID__k[ID_214,Intercept]"    "r_ID__k[ID_217,Intercept]"   
[265] "r_ID__k[ID_223,Intercept]"    "r_ID__k[ID_224,Intercept]"    "r_ID__k[ID_227,Intercept]"   
[268] "r_ID__k[ID_23,Intercept]"     "r_ID__k[ID_231,Intercept]"    "r_ID__k[ID_234,Intercept]"   
[271] "r_ID__k[ID_24,Intercept]"     "r_ID__k[ID_25,Intercept]"     "r_ID__k[ID_26,Intercept]"    
[274] "r_ID__k[ID_27,Intercept]"     "r_ID__k[ID_28,Intercept]"     "r_ID__k[ID_29,Intercept]"    
[277] "r_ID__k[ID_3,Intercept]"      "r_ID__k[ID_30,Intercept]"     "r_ID__k[ID_32,Intercept]"    
[280] "r_ID__k[ID_34,Intercept]"     "r_ID__k[ID_35,Intercept]"     "r_ID__k[ID_36,Intercept]"    
[283] "r_ID__k[ID_37,Intercept]"     "r_ID__k[ID_38,Intercept]"     "r_ID__k[ID_39,Intercept]"    
[286] "r_ID__k[ID_4,Intercept]"      "r_ID__k[ID_40,Intercept]"     "r_ID__k[ID_41,Intercept]"    
[289] "r_ID__k[ID_42,Intercept]"     "r_ID__k[ID_45,Intercept]"     "r_ID__k[ID_47,Intercept]"    
[292] "r_ID__k[ID_48,Intercept]"     "r_ID__k[ID_49,Intercept]"     "r_ID__k[ID_50,Intercept]"    
[295] "r_ID__k[ID_51,Intercept]"     "r_ID__k[ID_52,Intercept]"     "r_ID__k[ID_53,Intercept]"    
[298] "r_ID__k[ID_55,Intercept]"     "r_ID__k[ID_56,Intercept]"     "r_ID__k[ID_57,Intercept]"    
[301] "r_ID__k[ID_58,Intercept]"     "r_ID__k[ID_59,Intercept]"     "r_ID__k[ID_6,Intercept]"     
[304] "r_ID__k[ID_60,Intercept]"     "r_ID__k[ID_61,Intercept]"     "r_ID__k[ID_62,Intercept]"    
[307] "r_ID__k[ID_63,Intercept]"     "r_ID__k[ID_64,Intercept]"     "r_ID__k[ID_65,Intercept]"    
[310] "r_ID__k[ID_67,Intercept]"     "r_ID__k[ID_68,Intercept]"     "r_ID__k[ID_69,Intercept]"    
[313] "r_ID__k[ID_7,Intercept]"      "r_ID__k[ID_72,Intercept]"     "r_ID__k[ID_73,Intercept]"    
[316] "r_ID__k[ID_74,Intercept]"     "r_ID__k[ID_75,Intercept]"     "r_ID__k[ID_76,Intercept]"    
[319] "r_ID__k[ID_77,Intercept]"     "r_ID__k[ID_78,Intercept]"     "r_ID__k[ID_79,Intercept]"    
[322] "r_ID__k[ID_8,Intercept]"      "r_ID__k[ID_80,Intercept]"     "r_ID__k[ID_81,Intercept]"    
[325] "r_ID__k[ID_82,Intercept]"     "r_ID__k[ID_84,Intercept]"     "r_ID__k[ID_85,Intercept]"    
[328] "r_ID__k[ID_86,Intercept]"     "r_ID__k[ID_87,Intercept]"     "r_ID__k[ID_88,Intercept]"    
[331] "r_ID__k[ID_89,Intercept]"     "r_ID__k[ID_9,Intercept]"      "r_ID__k[ID_90,Intercept]"    
[334] "r_ID__k[ID_91,Intercept]"     "r_ID__k[ID_92,Intercept]"     "r_ID__k[ID_94,Intercept]"    
[337] "r_ID__k[ID_95,Intercept]"     "r_ID__k[ID_96,Intercept]"     "r_ID__k[ID_98,Intercept]"    
[340] "r_ID__k[ID_99,Intercept]"     "r_ID__I[ID_1,Intercept]"      "r_ID__I[ID_10,Intercept]"    
[343] "r_ID__I[ID_100,Intercept]"    "r_ID__I[ID_101,Intercept]"    "r_ID__I[ID_102,Intercept]"   
[346] "r_ID__I[ID_103,Intercept]"    "r_ID__I[ID_104,Intercept]"    "r_ID__I[ID_105,Intercept]"   
[349] "r_ID__I[ID_106,Intercept]"    "r_ID__I[ID_107,Intercept]"    "r_ID__I[ID_109,Intercept]"   
[352] "r_ID__I[ID_111,Intercept]"    "r_ID__I[ID_112,Intercept]"    "r_ID__I[ID_113,Intercept]"   
[355] "r_ID__I[ID_114,Intercept]"    "r_ID__I[ID_116,Intercept]"    "r_ID__I[ID_117,Intercept]"   
[358] "r_ID__I[ID_118,Intercept]"    "r_ID__I[ID_119,Intercept]"    "r_ID__I[ID_12,Intercept]"    
[361] "r_ID__I[ID_120,Intercept]"    "r_ID__I[ID_121,Intercept]"    "r_ID__I[ID_122,Intercept]"   
[364] "r_ID__I[ID_123,Intercept]"    "r_ID__I[ID_124,Intercept]"    "r_ID__I[ID_125,Intercept]"   
[367] "r_ID__I[ID_126,Intercept]"    "r_ID__I[ID_127,Intercept]"    "r_ID__I[ID_128,Intercept]"   
[370] "r_ID__I[ID_129,Intercept]"    "r_ID__I[ID_13,Intercept]"     "r_ID__I[ID_130,Intercept]"   
[373] "r_ID__I[ID_132,Intercept]"    "r_ID__I[ID_133,Intercept]"    "r_ID__I[ID_134,Intercept]"   
[376] "r_ID__I[ID_135,Intercept]"    "r_ID__I[ID_136,Intercept]"    "r_ID__I[ID_137,Intercept]"   
[379] "r_ID__I[ID_138,Intercept]"    "r_ID__I[ID_139,Intercept]"    "r_ID__I[ID_14,Intercept]"    
[382] "r_ID__I[ID_142,Intercept]"    "r_ID__I[ID_143,Intercept]"    "r_ID__I[ID_144,Intercept]"   
[385] "r_ID__I[ID_147,Intercept]"    "r_ID__I[ID_148,Intercept]"    "r_ID__I[ID_15,Intercept]"    
[388] "r_ID__I[ID_152,Intercept]"    "r_ID__I[ID_153,Intercept]"    "r_ID__I[ID_154,Intercept]"   
[391] "r_ID__I[ID_155,Intercept]"    "r_ID__I[ID_156,Intercept]"    "r_ID__I[ID_157,Intercept]"   
[394] "r_ID__I[ID_158,Intercept]"    "r_ID__I[ID_16,Intercept]"     "r_ID__I[ID_160,Intercept]"   
[397] "r_ID__I[ID_161,Intercept]"    "r_ID__I[ID_163,Intercept]"    "r_ID__I[ID_165,Intercept]"   
[400] "r_ID__I[ID_166,Intercept]"    "r_ID__I[ID_169,Intercept]"    "r_ID__I[ID_17,Intercept]"    
[403] "r_ID__I[ID_170,Intercept]"    "r_ID__I[ID_171,Intercept]"    "r_ID__I[ID_172,Intercept]"   
[406] "r_ID__I[ID_174,Intercept]"    "r_ID__I[ID_175,Intercept]"    "r_ID__I[ID_176,Intercept]"   
[409] "r_ID__I[ID_177,Intercept]"    "r_ID__I[ID_18,Intercept]"     "r_ID__I[ID_180,Intercept]"   
[412] "r_ID__I[ID_181,Intercept]"    "r_ID__I[ID_186,Intercept]"    "r_ID__I[ID_189,Intercept]"   
[415] "r_ID__I[ID_19,Intercept]"     "r_ID__I[ID_191,Intercept]"    "r_ID__I[ID_192,Intercept]"   
[418] "r_ID__I[ID_197,Intercept]"    "r_ID__I[ID_198,Intercept]"    "r_ID__I[ID_199,Intercept]"   
[421] "r_ID__I[ID_2,Intercept]"      "r_ID__I[ID_20,Intercept]"     "r_ID__I[ID_201,Intercept]"   
[424] "r_ID__I[ID_202,Intercept]"    "r_ID__I[ID_207,Intercept]"    "r_ID__I[ID_208,Intercept]"   
[427] "r_ID__I[ID_21,Intercept]"     "r_ID__I[ID_212,Intercept]"    "r_ID__I[ID_214,Intercept]"   
[430] "r_ID__I[ID_217,Intercept]"    "r_ID__I[ID_223,Intercept]"    "r_ID__I[ID_224,Intercept]"   
[433] "r_ID__I[ID_227,Intercept]"    "r_ID__I[ID_23,Intercept]"     "r_ID__I[ID_231,Intercept]"   
[436] "r_ID__I[ID_234,Intercept]"    "r_ID__I[ID_24,Intercept]"     "r_ID__I[ID_25,Intercept]"    
[439] "r_ID__I[ID_26,Intercept]"     "r_ID__I[ID_27,Intercept]"     "r_ID__I[ID_28,Intercept]"    
[442] "r_ID__I[ID_29,Intercept]"     "r_ID__I[ID_3,Intercept]"      "r_ID__I[ID_30,Intercept]"    
[445] "r_ID__I[ID_32,Intercept]"     "r_ID__I[ID_34,Intercept]"     "r_ID__I[ID_35,Intercept]"    
[448] "r_ID__I[ID_36,Intercept]"     "r_ID__I[ID_37,Intercept]"     "r_ID__I[ID_38,Intercept]"    
[451] "r_ID__I[ID_39,Intercept]"     "r_ID__I[ID_4,Intercept]"      "r_ID__I[ID_40,Intercept]"    
[454] "r_ID__I[ID_41,Intercept]"     "r_ID__I[ID_42,Intercept]"     "r_ID__I[ID_45,Intercept]"    
[457] "r_ID__I[ID_47,Intercept]"     "r_ID__I[ID_48,Intercept]"     "r_ID__I[ID_49,Intercept]"    
[460] "r_ID__I[ID_50,Intercept]"     "r_ID__I[ID_51,Intercept]"     "r_ID__I[ID_52,Intercept]"    
[463] "r_ID__I[ID_53,Intercept]"     "r_ID__I[ID_55,Intercept]"     "r_ID__I[ID_56,Intercept]"    
[466] "r_ID__I[ID_57,Intercept]"     "r_ID__I[ID_58,Intercept]"     "r_ID__I[ID_59,Intercept]"    
[469] "r_ID__I[ID_6,Intercept]"      "r_ID__I[ID_60,Intercept]"     "r_ID__I[ID_61,Intercept]"    
[472] "r_ID__I[ID_62,Intercept]"     "r_ID__I[ID_63,Intercept]"     "r_ID__I[ID_64,Intercept]"    
[475] "r_ID__I[ID_65,Intercept]"     "r_ID__I[ID_67,Intercept]"     "r_ID__I[ID_68,Intercept]"    
[478] "r_ID__I[ID_69,Intercept]"     "r_ID__I[ID_7,Intercept]"      "r_ID__I[ID_72,Intercept]"    
[481] "r_ID__I[ID_73,Intercept]"     "r_ID__I[ID_74,Intercept]"     "r_ID__I[ID_75,Intercept]"    
[484] "r_ID__I[ID_76,Intercept]"     "r_ID__I[ID_77,Intercept]"     "r_ID__I[ID_78,Intercept]"    
[487] "r_ID__I[ID_79,Intercept]"     "r_ID__I[ID_8,Intercept]"      "r_ID__I[ID_80,Intercept]"    
[490] "r_ID__I[ID_81,Intercept]"     "r_ID__I[ID_82,Intercept]"     "r_ID__I[ID_84,Intercept]"    
[493] "r_ID__I[ID_85,Intercept]"     "r_ID__I[ID_86,Intercept]"     "r_ID__I[ID_87,Intercept]"    
[496] "r_ID__I[ID_88,Intercept]"     "r_ID__I[ID_89,Intercept]"     "r_ID__I[ID_9,Intercept]"     
[499] "r_ID__I[ID_90,Intercept]"     "r_ID__I[ID_91,Intercept]"     "r_ID__I[ID_92,Intercept]"    
[502] "r_ID__I[ID_94,Intercept]"     "r_ID__I[ID_95,Intercept]"     "r_ID__I[ID_96,Intercept]"    
[505] "r_ID__I[ID_98,Intercept]"     "r_ID__I[ID_99,Intercept]"     "lp__"                        

lobe.n.best.fitted.summary <- summary(lobe.n.best.fitted)$summary
write.csv(lobe.n.best.fitted.summary, file = "/Users/seungmokim/414_Growth_Model/summary/growth_model_phenotypes/lobe_number/lobe.n.best.fitted.summary.csv")

r_ID__Hmax

datalist0 = list()
for (samples in sample.list) {
  ttt0 <- paste0("r_ID__Hmax[",samples,",Intercept]")
  dat0 <- extract(lobe.n.best.fitted, ttt0, permuted=TRUE)
  datalist0[[ttt0]] <- dat0
}
list.data0 = do.call(rbind, datalist0)
Hmax.ID <- data.frame(matrix(unlist(list.data0), nrow = 166, byrow = T))
row.names(Hmax.ID) <- paste0("r_ID__Hmax[",sample.list,",Intercept]")
head(Hmax.ID)

r_ID__k

datalist1 = list()
for (samples in sample.list) {
  ttt1 <- paste0("r_ID__k[",samples,",Intercept]")
  dat1 <- extract(lobe.n.best.fitted, ttt1, permuted=TRUE)
  datalist1[[ttt1]] <- dat1
}
list.data1 = do.call(rbind, datalist1)
k.ID <- data.frame(matrix(unlist(list.data1), nrow = 166, byrow = T))
row.names(k.ID) <- paste0("r_ID__k[",sample.list,",Intercept]")
head(k.ID)

r_ID__I

datalist2 = list()
for (samples in sample.list) {
  ttt2 <- paste0("r_ID__I[",samples,",Intercept]")
  dat2 <- extract(lobe.n.best.fitted, ttt2, permuted=TRUE)
  datalist2[[ttt2]] <- dat2
}
list.data2 = do.call(rbind, datalist2)
I.ID <- data.frame(matrix(unlist(list.data2), nrow = 166, byrow = T))
row.names(I.ID) <- paste0("r_ID__I[",sample.list,",Intercept]")
head(I.ID)

estimates of parameters

par.list = c("b_Hmax_Intercept", "b_k_Intercept", "b_I_Intercept", "b_Hmin_Intercept", "sd_ID__Hmax_Intercept",    "sd_ID__k_Intercept","sd_ID__I_Intercept","sigma", "lp__")
datalist3 = list()
for (pars in par.list) {
  dat3 <- extract(lobe.n.best.fitted, pars, permuted=TRUE)
  datalist3[[pars]] <- dat3
}
list.data3 = do.call(rbind, datalist3)
pars.est <- data.frame(matrix(unlist(list.data3), nrow = 9, byrow = T))
row.names(pars.est) <- par.list
pars.est

combine estimates of all parameters

pars.estimate.lobe.n.fitted <- rbind(pars.est, Hmax.ID, k.ID, I.ID)
head(pars.estimate.lobe.n.fitted, 10)

write.csv(pars.estimate.lobe.n.fitted, file="pars.estimate.lobe.n.fitted.csv")